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Abstract 

We formulate a lattice Boltzmann algorithm which solves the hydrodynamic equations of motion for nematic 
liquid crystals. The applicability of the approach is demonstrated by presenting results for two liquid crystal 
devices where flow has an important role to play in the switching. 
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1. Introduction 

In this paper we describe a lattice Boltzmann al- 
gorithm which solves the hydrodynamic equations 
of motion for nematic liquid crystals. Liquid crys- 
tals are widely used in display devices and, as ex- 
amples of applications of the algorithm, we con- 
sider two model devices where flow has an impor- 
tant role to play in the switching. 

Liquid crystals are fluids made up of rod- 
shaped molecules[l]. At high temperatures and 
low concentrations the angular distribution of the 
molecules is isotropic. However, as the tempera- 
ture is lowered or the concentration increased, a 
liquid crystal can undergo a phase transition to a 
state where the molecules tend to align parallel, 
the so-called nematic phase. Liquid crystals ex- 
hibit complicated, non-Newtonian flow behaviour 
because of the coupling between this molecular 
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structure and the flow field [2]. Examples include 
shear thinning and non-equilibrium phase transi- 
tions such as banding under shear [3]. 

Another consequence of the ordering of nematic 
liquid crystals is their widespread use in optical 
display devices such as those on digital watches or 
calculators. In a typical display device, the liquid 
crystal is confined between two plates a few mi- 
crons apart. The molecular configuration on the 
plates is fixed. When an electric field is switched 
on molecules in the bulk align in the direction 
preferred by the field. After switching off the field, 
long-range elastic interactions ensure that the 
molecules reorient themselves in the direction pre- 
ferred by the surfaces. These devices can be used 
as displays because different liquid crystal orien- 
tations have different optical properties. Flow can 
be important in device operation because it can 
control the speed of switching or even select the 
states between which switching can occur. 

Modelling flow in non-Newtonian fluids such as 
liquid crystals is difficult because of the need to 
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incorporate the coupling between the microscopic 
structure and the flow. The hydrodynamic equa- 
tions of motion can be extremely complex and nu- 
merical solutions demanding. A method which has 
recently proved very successful in modelling com- 
plex fluids is the lattice Boltzmann approach[4]. 
Lattice Boltzmann simulations solve the hydro- 
dynamic equations of motion while inputting suf- 
ficient, albeit generic, molecular information to 
model the important physics of a given fluid. This 
is often done by imposing a Landau free energy 
functional, so the fluid evolves to a known ther- 
modynamic equilibrium [5]. The method has been 
successfully used to investigate domain growth in 
binary mixtures and liquid-gas systems [6] , the flow 
of binary mixtures in porous media, and ordering 
in amphiphilic fluids [7]. 

Here we construct a lattice Boltzmann algorithm 
to solve the hydrodynamic equations of motion for 
nematic liquid crystals. We follow Beris and Ed- 
wards [8] and consider a rather general formalism 
of the hydrodynamics written in terms of a ten- 
sor order parameter, Q. This approach lends itself 
to a lattice Boltzmann interpretation: the partial 
density distribution functions which are the usual 
variables in the simulations are supplemented by 
a set of tensors related to Q. Backflow, the hydro- 
dynamics of topological defects and the possibility 
of transitions between the nematic and isotropic 
phases appear naturally within the formalism. 

The next section of the paper summarises the 
Beris-Edwards equations of motion for liquid crys- 
tal hydrodynamics. Sec. 3 describes the lattice 
Boltzmann algorithm. General viscosity terms, 
needed to match to experimental viscosities, are 
included. The approach is applied to switching in 
a Frederiks cell in Sec. 4 and, in Sec. 5, used to 
confirm that backflow can influence state selection 
in a bistable liquid crystal device. 



2. Hydrodynamic equations of motion 

Ordering in a nematic liquid crystal can be de- 
scribed in terms of a tensor order parameter Q 
which is symmetric and traceless. Expanding the 
free energy in terms of the order parameter gives 



the Landau-de-Gennes expression[l] 
$= f dV {F b - F E + F d } + f dS{F a }, (1) 

JV JdV 

where 
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(Greek subscripts represent Cartesian directions 
and the usual summation over repeated indices is 
assumed.) The bulk free energy terms Ft describe a 
liquid crystal with a first-order, isotropic-nematic 
transition at 7 = 2. 7 [9]. Contributions to the free 
energy which arise from an imposed electric field 
E arc included in F E ■ Fd describes the elastic free 
energy [8] and F a is a surface free energy which fixes 
a preferred orientation Q° for the surface director 
field. 

Q evolves according to a convection-diffusion 
equation [8] 

(d t +u - V)Q - S(W, Q) = FH (3) 

where u is the bulk fluid velocity and V is a col- 
lective rotational diffusion constant. The term on 
the right-hand side of Eqn. (3) describes the relax- 
ation of the order parameter towards the minimum 
of the free energy 

H = -| + (I/3)T r {|}. ,4) 

The order parameter distribution can be both ro- 
tated and stretched by flow gradients. This is de- 
scribed by the term on the left-hand side 

S(W,Q) = (£D + ft)(Q + I/3) + (5) 
(Q + 1/3) (^D - ft) - 2£(Q + I/3)Tr(QW) 

where D = (W + W T )/2 and ft = (W - W T )/2 
are the symmetric part and the anti-symmetric 
part respectively of the velocity gradient tensor 
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W a p = dpu a and £ is related to the aspect ratio of 
the molecules. 

The flow of the liquid crystal fluid obeys the 
continuity and Navier-Stokes equations 



dtp + d a pu a = 0, 



pd t U a + pilpdpUa = 8p(T a p + 
pTf, 



(6) 



(7) 
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where 
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(8) 

is the stress due to the nematic order [8] and the 
external electric field. The latter is represented by 
the <JM,ap Maxwell stress tensor [10]. 



where the collision operators are taken to have the 
form of a single relaxation time Boltzmann equa- 
tion, together with a forcing term 



C fi {x,t,{fi}) = 

fi(x,t)-f?(x,t,{fi}) 



(11) 



+ Pi(x,t,{fi}), 



C Gi (x,t,{G i }) = 

G l (x,t)-G^(x,t,{G t }) 



+ hi(x,t,{Gi}). 



f* and G* in equations (10) and (11) are first or- 
der approximations to fc(x + etAt,t + At) and 
Gi(x + e,At, t + At). They are introduced to re- 
move lattice viscosity terms to second order and 
they give improved stability. 

The form of the equations of motion and ther- 
modynamic equilibrium follow from the choice of 
the moments of the equilibrium distributions ff 
and G° and the driving terms pi and hi . ff is con- 
strained by 



i i 
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(12) 



3. Lattice Boltzmann algorithm 

A lattice Boltzmann scheme which reproduces 
Eqns. (3), (6) and (7) to second order can be de- 
fined in terms of two distribution functions fi(x) 
and Gi(x) where i labels lattice directions from 
site x. Physical variables are related to the distri- 
bution functions by 

/° = ^/*> pu a = Y-f ieia > Q = X] Gl 

i i i 

(9) 

The distribution functions evolve in a time step 
At according to 



(fi(x + e t At, t + At)- fi(x, t))/At : 



(10) 



- [C fi (x, t, {fi}) + C fi (x + ei At, t + At, {/*})] : 

(Gi(x + e t At, t + At)- Gi{x, t))/At = 
C Gl (x, t, {GO) + C Gl (x + e t At, t + At, {G*}) 



where the zeroth and first moments are chosen to 
impose conservation of mass and momentum. The 
second moment of / controls the symmetric part 
of the stress tensor, whereas the moments of Pi 



i i 



(13) 



impose the antisymmetric part of the stress tensor. 

For the equilibrium of the order parameter dis- 
tribution we choose 



^G° = Q, ^G°e ia = Q« c 

i i 

^ G°e ia eip = Qu a up. 



(14) 



This ensures that the fluid minimizes its free energy 
at equilibrium and that it is convected with the 
flow. Finally the evolution of the order parameter 
is most conveniently modelled by choosing 
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Fig. 2. Angle (9, with the horizontal, of the director at the 
centre of the cell as a function of time. The system starts 
in the fully relaxed, zero voltage state and the voltage 
is turned on immediately at the far left of the plot. The 
voltage is turned off at t = 0. 



Fig. 1. The two stable states of a Frederiks cell. At zero 
voltage the liquid crystal prefers the orientation imposed 
by the boundaries. When the electric field is turned on the 
nematic prefers to align along the field direction. 

]T^=rH(Q) + S(W,Q), 

i 

^ he ia = (y] hj)u a . (15) 

i i 

Conditions (12)-(15) can be satisfied by taking /?, 
G", hi, and pi as polynomial expansions in the ve- 
locity as is usual in lattice Boltzmann schemes [4]. 
A second order Chapman-Enskog expansion for 
the evolution equations (10) incorporating the con- 
ditions (12)-(15) leads to the equations of mo- 
tional]. 



4. Frederiks cell 

One of the simplest liquid crystal devices is the 
Frederiks cell depicted in Fig. 1. When the volt- 
age is increased above a critical value V c switching 
from the zero voltage state occurs rather rapidly. 
The factor limiting the speed at which the device 
can be updated is the switching from the voltage 
aligned state to that preferred by the boundaries. 
This switching involves both relaxation of the free 
energy and viscous dissipation associated with hy- 
drodynamic modes[l]. Our aim is to quantitatively 
model this phenomenon for a given set of elastic 
constants and viscosities [12]. 



Figure 2 shows 9, the director angle [13] mea- 
sured in the centre of the channel, as indicated 
in Fig. 1, as a function of time as the voltage is 
switched on and then off. The parameters used 
for the simulation are listed in footnote [14] and 
were chosen to correspond to commercially rele- 
vant materials. For comparison, we show relax- 
ation data, provided by Sharp Laboratories of Eu- 
rope, taken from a simulation using a commercial 
code which solves the Ericksen-Leslie[15] equations 
in one dimension. (For a similar, but not iden- 
tical, Ericksen-Leslie calculation and comparison 
to experimental data see Ref. [16].) For the most 
part, there is very good agreement between the two 
methods. 

The "blip" seen in the plot when the voltage 
is turned off is a real effect known as the optical 
bounce, due to its optical signature. After the re- 
moval of the electric field, the changing director 
field induces a hydrodynamic flow. This, in turn, 
couples back to the director field and causes it 
to momentarily move in the "wrong" direction. 
Viscous dissipation in the fluid rapidly damps out 
the flow and the director can then relax toward 
the minimum of the free energy. The bounce we 
see is somewhat higher than that observed in 
the Ericksen-Leslie simulations. However recent 
attempts to exactly match the quantitative fea- 
tures of the optical bounce to solutions of the 
Ericksen-Leslie equations have suggested that 
they always underestimate its value relative to 
experiments [16]. Thus, the additional details in 
the tensor model may be better reproducing the 
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Fig. 3. The two stable states (a) horizontal and (b) twist 
of the surface switching device considered in Sec. 5. In (b) 
the molecular orientation twists to perpendicular to the 
plane of the figure at the middle of the sample. 

experimental results. 



5. Flow induced surface switching 

A current aim of the device industry is to con- 
struct bistable displays which can remain in two 
different (meta)stable states in zero field. This 
would, for example, lead to significant power sav- 
ing in infrequently updated displays or could be 
used in the construction of smart cards which rely 
on an external power source. 

In this paper we will consider one possible mech- 
anism for accessing two zero-field states, flow in- 
duced surface switching [17] . The two states are the 
horizontal and the twist states shown in Fig. 3. 
The liquid crystal is confined between two plates. 
The top (bottom) plate has strong (weak) anchor- 
ing. The horizontal state can be switched to the 
vertical state by an electric field. Then, if the field 
is abruptly switched off, one might expect the sys- 
tem to relax back to the horizontal state. The re- 
laxation commences at the boundaries, where the 
molecules start to rotate in the same direction at 
the top and bottom surfaces. If one does not con- 
sider the effect of the back-flow (equivalent to us- 
ing a Ginzburg-Landau equation to describe the 
dynamics), then this is indeed the shown 
by the simulation results in Fig. 4. (The simula- 
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Fig. 4. Snapshots of the director configuration of the surface 
switching device after the electric field is switched off, but 
back-flow effects are not incorporated in the model. The 
(a) vertical state relaxes through (b) an intermediate state 
to (c) the horizontal state. The snapshots correspond to 0, 
7 and 15 ms. At the top surface the molecules are strongly 
pinned at a 20° pretilt to the horizontal axis. Pinning at 
the bottom surface is weak with a 0° pretilt. 



tion parameters are summarized in [18]. Back-flow 
is switched off by setting the velocity to zero in 
Eqn. (3). 

If, however, back-flow is included, a horizontal 
flow is created starting from the top plate and grad- 
ually reaching the bottom surface. This horizontal 
flow causes the molecules to rotate in the opposite 
direction at the bottom surface. As a result the ver- 
tical state transforms to a bend state, as shown in 
Fig. 5. There is an energy argument which suggests 
that the bend state is unstable and quickly relaxes 
to a twist state [17] if the Frank elastic constant 
if 22 is sufficiently small compared to K\\ and K^. 

Hence in a bistable device, if the field is switched 
off gradually, the effect of back-flow is weak, and 
the liquid crystal settles in the horizontal state. 
If the field is switched off abruptly, then it settles 
in the twist state. When a field is reapplied both 
zero field configurations switch back to the vertical 
state. 



5 



///// / / / / / s s s y s 

I I I I I / / / / / 

I I I I I ///// / / / s / 

\ \ \ \ \ ///// / / / / / 

I I I I I ///// ///// 

I I I I I ///// ///// 

I I I I I ///// ///// 

I I I I I I I I I I ///// 

I I I I I I I I I I I I I I I 

I I I I I I I I I I i i I I I 

I I I I I I I I I I I I 1 I I 

I I I I I I I I I I I I I I I 

I I I I I I I I I I i I I I ! 

I I I I I I I I I I I I I I I 

I I I I I I I 1 I I I I I I I 

I I I I I I I 1 I I I 1 I 1 I 

I I I I I I I I I I \ 

I I I I I I I I I I \ \ \ \ \ 

I I I I I 11111 \ \ \ \ \ 

I I I I I 11111 \ \ \ \ \ 

III 11111 \\\\\ 

III \ \ \ \ \ \ \ \ \ \ 

III \ \ \ \ \ \\\\\ 

III \\\\\ \\\\\ 

I I I I I \\\\\ \\\\\ 

III \\\\\ \ \ N \ N 



(a) (b) (c) 

Fig. 5. Snapshots of the director configuration of the sur- 
face switching device after the electric field is switched off, 
and back-flow effects arc incorporated in the model. The 
(a) vertical state relaxes through (b) an intermediate state 
to (c) the bend state. The snapshots correspond to 0, f 
and 6 ms. 

6. Conclusion 

In conclusion, we have derived and implemented 
a lattice Boltzmann algorithm for liquid crystal 
hydrodynamics. This opens the way to investigate 
a wide range of physical phenomena which result 
from the coupling between the director field and 
the flow, many of which have been examined only 
indirectly or with severe approximations in the 
past. Two examples have been presented showing 
the relevance of flow to switching in liquid crystal 
devices. 
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